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ABSTRACT 
Methods  are  developed,  and  Fortran  63  CODAP  computer  programs 
are  demonstrated,  to  generate  random  numbers  from  the  uniform, 
normal  (including  multivariate  normal),  Poisson,  and  exponential 
probability  distributions.  Various  statistical  tests  are  described 
and  the  results  of  the  application  of  these  tests  to  the  generators 
are  tabulated.  A  general  method  for  generating  random  numbers  from 
a  large  class  of  distributions  is  described.  The  methods  of 
generation  are  optimized  to  provide  an  accurate  generator  while 
producing  numbers  at  a  maximum  rate.  The  uniform  generator  that  is 
used  as  a  basis  for  the  other  generators  is  of  the  congruential 
type  and  is  capable  of  generating  1800  numbers  per  second. 
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1.   Introduction. 

The  need  for  a  rapid  reliable  source  of  random  numbers  from 
a  prescribed  distribution  is  well  recognized,  and  is  likely  to 
become  even  more  pressing  in  military  circles  in  view  of  the 
Department  of  Defense's  increased  interest  in  Operations  Analysis. 
This  thesis  presents  methods  and  programs,  some  old,  some  new, 
for  generating  random  numbers  from  various  specified  distributions. 
Some  statistical  tests  of  the  programs  and  their  results  are 
described. 

All  methods  assume  a  uniform  random  number  generator  is 
available.  A  thesis  by  Barron  [lj  has  a  good  bibliography  on 
this  subject.  The  generator  used  here,  and  described  more  completely 
in  Section  3>  is  of  the  mixed  congruential  type.  While  some  uniform 
generators  may  have  advantages  over  the  one  used,  this  one  seems 
to  perform  very  well,  at  the  same  time  as  being  as  fast  as  any 
demonstrated.  Since  this  generator  is  used  as  a  basis  for  all  the 
others  it  should  be  remembered  that  no  generator  can  be  considered 
•perfect1,  especially  in  the  continuous  distribution  case,  since  the 
computer  is  limited  to  a  finite  set  of  possible  numbers.  However, 
for  practical  purposes  this  inaccuracy  is  not  important.  Other 
sources  of  inaccuracy  c an  be  important,  however.  The  numbers 
generated  must  have  the  property  of  randomness,  and  must  faithfully 
represent  the  desired  distribution.  These  properties  are  measured 
by  testing  samples  of  the  generated  numbers.  Some  of  the  methods  of 
generating  numbers  in  this  paper  theoretically  provide  an  exact 


transformation  from  the  uniform  distribution  to  the  desired 

distribution.  An  example  of  this  is  the  half -Gaussian  method 

1 
of  generating  normal  random  numbers  .  If  we  assume  that  the 

uniform  numbers  are  accurate  then  we  are  led  to  the  conclusion 

that  the  normal  numbers  are  also  accurate  and  that  tests  of  these 

numbers  would  be  superfluous.  However,  tests  are  performed  to 

assure  that  the  uniform  numbers  were  so  good  that  they  did  not 

bias  the  derived  distribution.  Other  techniques  such  as  Marsaglia's 

technique  (see  Section  h)   for  generating  normal  random  numbers 

are  in  a  sense  curve-fitting  techniques  and  only  provide  a 

controllably  good  approximation  to  the  real  distribution.  The 

advantage  of  the  approximation  techniques  is  in  the  far  greater 

speed  with  which  they  may  provide  the  desired  numbers.  The  user 

who  needs  to  draw  numbers  from  a  distribution  he  suspects  is 

normal,  with  inaccurately  measured  mean  and  hypothesized  variance, 

is  not  in  need  of  a  generator  that  is  accurate  in  the  sixth  decimal 

place,  however,  he  still  would  like  to  be  assured  that  having 

assumed  a  distribution  and  its  parameters  he  will  be  able  to  generate 

numbers  with  the  appropriate  shape  and  with  good  properties  of 

randomness. 


xhe  phrase  '  normal  random  numbers  ■  and  other  like  phrases 
should  be  read  as  'numbers  distributed  as  if  coming  from  the  normal 
distribution. ' 


2.   Testing  the  Generators. 
2.1  General  Discussion. 

The  uniform  generator  chosen  here  has  been  tested  extensively 
by  others.  The  tests  that  have  been  applied  include  frequency 
tests,  serial  tests,  moment  tests,  poker  tests,  gap  tests,  and 
many  others.  As  mentioned  in  the  introduction,  this  generator  is 
as  good  as  can  be  found,  considering  the  requirement  for  speed. 
However,  the  derived  distributions  will  be  tested  to  overcome  any 
doubts  there  may  be  about  the  transformation.  The  numbers  will 
be  tested  mainly  to  measure  the  faithfulness  with  which  they 
represent  the  derived  distribution.  The  randomness  is  provided 
by  the  uniform  generator.  If  the  randomness  is  not  satisfactory 
then  the  uniform  generator  must  be  blamed,  not  the  transformation. 
A  slower  but  more  satisfactory  method  is  available  if  the  need  is 
felt.    Thus  the  tests  used  here,  the  moment  test,  the  hypothesis 
tests  on  the  mean  and  variance,  and  the  Kolmogorov-Smirnov  goodness 
of  fit  test,  are  not  designed  to  detect  special  types  of  non- 
randomness,  such  as  is  detected  by  the  poker  test  and  other  similar 
tests. 

Martin  Greeriberger  has  written  an  interesting  article  [17]  on 
this  subject.  He  presents  the  results  of  an  investigation  by 
Joseph  Lach  [l6Q  at  Yale  University  in  which  Lach  showed  that  the 
congruential  method  of  generating  uniform  random  numbers  has  a 


"See  page  22 » 


predictable  non-randomness  in  second  order  serial  correlation,, 
The  lesson  is  clear.  Having  found  an  undersirable  feature  in  a 
generator,  it  is  generally  possible  to  modify  the  generator  to 
eliminate  the  feature,  however,  we  can  be  sure  that  we  have 
introduced  another  aberration  of  some  kind,  even  though  its  form 
may  be  hard  to  determine.  When  the  user  asks,  "Is  this  generator 
good  enough?",  the  obvious  retort  is  "good  enough  for  what?" 
No  one  generator  is  suitable  for  all  applications,  but  the  generator 
used  here  will  be  good  enough  for  most.  If  the  user  thinks  that 
this  is  not  so  in  his  application,  he  has  at  his  resources  the 
modifying  methods  of  Marsaglia  Id}   or  Lach  with  the  penalty  of 
longer  generation  times. 
2.2  Moment  tests 

The  first  four  moments  are  calculated  and  are  compared  with 
the  theoretical  moments  for  each  of  the  distributions.  A  more 
appealing  statistic  than  the  sample  moment  may  be  the  unbiased 
estimator  of  the  moment,  however  for  large  sample  sizes  such  as 
are  used  here,  this  varies  very  little  from  the  sample  moment,  and 
the  sample  moment  is  easier  to  handle  in  other  uses-  such  as  hypothesis 
testing.  The  unbiased  estimator  of  the  variance  is: 

The  statistic  used  here  is  the  sample  second  moment: 

All  numerations  iro  on  the  index  '4.'  which  runs  from  1  to  N, 
the  sample  size,  unless  otherwise  indicated. 


The  first  moment  is  the  mean: 
The  third  and  fourth  moments  aret 

2.3    Hypothesis  tests  on  the  mean  and  variance. 

The  availability  of  large  sample  sizes  is  used  in 
designing  tests  on  the  mean  and  variance.  The  cental  limit 
theorem  is  used  where  possible  to  simplify  the  test  procedures. 

2.3.1  Test  on  the  mean* 

This  test  is  applied  to  the  normal  distribution.  The 
hypothesis  is  that  the  mean  is  zero;;  the  alternate  hypothesis 
is  that  the  mean  is  not  zero#  The  test  is  performed  by 
calculating  the  statistic  Y,  where 

Y=  In* 

The  decision  rule  at  the  alpha  level  becomes: 

Accept  the  hypothesis  when   -  K°V^  <•  I  <  KK/L 

At  the  10$  level  this  rule  becomes: 

Accept  the  hypothesis  when   —|ifrS<  Y<   LfrS 

The  justifications  for  using  this  test  in  preference  to  the  fT*  test 

are  that  the  variance  can  be  assumed  to  be  one,  and  that  a  large 

sample  size  is  being  used. 

2.3.2  Test  on  the  variance. 

Again  for  the  normal  distribution,  the  hypothesis  is  that 


the  variance  is  onej  the  alternate  hypothesis  is  that  the 
variance  is  not  one.  The  test  is  performed  by  calculating  the 
statistic  Z,  where     ^  -   $—  (.N-3JO 

The  decision  rule  at  the  10$  level  becomes s 

Accept  the  hypothesis  when   —  l.t^S  <  "2:  <  Lt^-'S. 
2.U    The  Kolmogorov-Smirnov  goodness  of  fit  test. 

In  Massey's  discussion  of  this  test  [6]    ,  he  presents  evidence 
to  indicate  that  this  test  may  in  many  circumstances  be  better 
than  the  more  usual  chi -squared  goodness  of  fit  test.  To  test 
the  hypothesis  that  a  sample  of  size  N  comes  from  theorized 
distribution,  the  cumulative  step  function  SN(x)  is  formed. 

SN(x)=k/n 
where  k  is  the  number  of  observations  less  than  or  equal  to  x. 
The  selection  of  x  is  arbitrary  within  certain  limits.  In  this 
paper  x  is  chosen  so  that  there  are  either  twenty  or  fifty  equal 
intervals  spanning  the  sample  space,  SN(x)  is  compared  with  the 
theoretical  value  of  the  cumulative  distribution,  F(x),  The 
maximum  difference  d  is  calculated, 

d=max|F(x)-SN(x)| 
Tables  due  to  Smirnov  [7  J  give  certain  critical  points  of  the 
distribution  for  various  sample  sizes.  For  sample  sizes  over 
35,  and  at  the  10$  level  of  significance,  if 

d/n<1.22/JT 


the  sampled  distribution  is  accepted  as  the  hypothesized 

distribution* 

2.5    The  chi -squared  goodness  of  fit  test. 

This  test  was  developed  by  Pearson  and  represents  the 
earliest  non-parametric  decision-making  test  in  statistics. 
Partly  because  of  the  length  of  time  the  test  has  been  in 
existence,  many  drawbacks  to  the  test  have  been  noted,  however 
it  still  stands  as  a  useful  and  much  used  test.  For  this  test 
the  sample  space  has  been  divided  into  k  intervals  and  the  number 
of  sample  observations  in  each  interval  is  noted. 

x    _    bi"^i 

where  n\  is  the  observed  number  of  sample  observations  in  each 
interval,  and  m-  is  the  expected  number.  Pearson  showed  that 
X  has  the  chi -squared  distribution  with  k-1  degrees  of  freedom. 
The  decision  rule  at  the  alpha  level  becomes:   if  X  —  -^\_|  (*' 
accept  the  hypothesis  that  the  distribtuion  is  as  postulated.  There 
are  several  problems  associated  with  the  application  of  this  test. 
How  should  the  interval  size  be  determined?  How  many  intervals 
should  there  be?  Mann  and  Wald  [Hi]  studied  this  problem  and 
formulated  a  criterion  for  the  selection  of  k,  Williams  [l5j  notes 
that  this  criterion  is  not  particularly  sensitive  to  even  a  reduction 
by  a  factor  of  two  in  the  number  of  intervals.  The  use  of  equal 
probability  intervals  vice  equal  length  intervals  is  also  recommended. 
However,  the  basis  for  this  recommendation  is  unclear.  It  is 
agreed  that  very  low  probability  intervals  such  as  would  occur 


in  the  tails  in  an  equal  interval  length  division  of  the  normal 
density  function  should  be  avoided.  The  test  is  applied  here 
using  equal  length  intervals.  Low  probability  intervals  are 
avoided  by  'pooling'  several  intervals  until  the  probability  is 
of  the  same  order  of  magnitude  as  in  the  other  intervals* 
2.6  Scatter  diagrams 

The  best  type  of  test  to  apply  to  the  generator  initially 
is  some  type  of  scatter  diagram.  The  scatter  diagram  can  often 
immediately  give  an  intuitive  idea  as  to  whether  the  generator  is 
behaving  properly.  In  fact  the  scatter  diagram  can  be  a  very 
powerful  tool  for  rejecting  a  generator-  more  sophisticated 
techniques  are  needed  to  accept  the  generator,  however.  The 
scatter  diagram  that  was  used  here  was  constructed  by  plotting 
the  first  number  generated  versus  the  second,  the  third  versus 
the  fourth,  and  so  on.  This  type  of  plot  will  also  enable  us  to 
look  for  correlations  similar  to  those  found  by  Lach  [l8J  and 
discussed  further  in  Section  2.1. 


3.     The  Uniform  Distribution, 
3.1    Disbribution  characteristics. 

It  is  desired  to  generate  numbers  such  that: 

The  first  four  moments  ares  Ml  =  1/2  :  M2  =  lA2  ',   M3  =  0  j 
MU  =  1/80.  The  uniform  distribution  is  usually  specified  in 
terms  of  its  interval  -  uniform  on  the  interval  from  z  ero  to  one 
(  denoted  here  by  U  (0,1)).  In  the  general  case  U  (A,B),  a  simple 
transformation  from  U  (0,1)  iss 

URAB  =  (UR01)(B-A)  +  A 
where  UR01  is  the  number  provided  by  the  U  (0,1)  generator, 
and  URAB  is  the  random  number  uniform  on  the  interval  (A,B). 
3 #2    Methods  of  generation. 

tt 

3.2.1  Many  techniques  have  been  used  over  the  years.  For  some 
particular  applications  such  a  method  as  table  look-up  may  be 
suitable.  However  for  our  purposes  what  is  desired  is  a  rapid, 
'accurate'  method  for  the  computer  to  produce  a  practically 
inexhaustible  supply  of  numbers. 

3.2.2  An  early  computational  scheme  was  called  the  mid-square 
method.  In  this  procedure  two  starting  values,  say  Al  and  A2, 
are  multiplied  together;  the  middle  set  of  bits  (usually  2U)  are 
extracted  as  the  third  random  number  A3$  then  A2  and  A3  are 
multiplied  together  and  the  algorithm  is  continued  in  a  similar 
fashion.  This  method  has  performed  well  in  many  tests  but 


unfortunately  degenerates  to  all  zeros  in  a  relatively  short 

time. 

3.2.3  An  improved  computational  scheme  was  tested  extensively 

by  Hull  and  Dobell  [3]  .  This  method  forms  a  series  of  numbers 

Ai,  where 

A7 

AR+1  =  BAR  +  ^      modulus  2 
B  and  C  are  constants  to  be  selected. 

3.2 .U  A  special  form  of  this  generator  where  C  =  0  is  the 
generally  recommended  form.  C  is  set  equal  to  zero  because  it 
does  not  improve  the  characteristics  of  the  generator  and  it  adds 
to  computation  time.  The  selection  of  B  and  the  starting  number  X 
has  been  the  subject  of  considerable  study.  Barron  and  others 
have  shown  that: 

XR+1  =  t2-19   +  3)XR  modulus  2^ 
where  X  is  either  1,  or  2^  -1,  or  any  number  naturally  generated  in 

the  sequence,  is  an  excellent  generator.  Since  this  generator  had 

been  tested  extensively  previously  no  attempt  was  made  to  test  it 

rigorsly  although  several  interesting  characteristics  were  noted. 

Some  of  these  are  included  here.  The  generator  was  run  down  through 

7     fi 
the  first  10  to  10  numbers.  A  number  at  the  end  of  this  sequence 

was  extracted  for  possible  alternate  use  as  a  starting  number.  It 

can  be  found  in  Appendix  I.  A  graph  of  the  mean  of  the  first  10000 

times  i  numbers ,  where  i  runs  from  1  to  100  is  plotted  against  the 

index  i.  This  plot  is  compared  with  curves  of  ■  K  qc/  lOOOOi  +  0.^0000 

versus  i.  The  more  our  data  plots  between  these  curves  the  better. 

10 


We  would  expect  it  to  be  between  the  curves  90$  of  the  time. 
As  the  following  graph  shows  our  generator  does  not  quite 
live  up  to  this  expectation.  A  scatter  diagram  of  the  type 
suggested  by  Lach  [l8]  was  plotted.  The  graph  on  page  13 
consists  of  3000  points  constructed  as  described  in  Section  2  .6. 

3.2.5  George  Marsaglia  and  M.  Donald  MacLaren  [8]  at  Boeing 
Scientific  Research  Laboratories  suggest  that  the  combination  of 
two  generators  will  produce  a  superior  random  number  of  generator. 
They  have  tested: 

AR+1  =  (217  +  3)AR  modulus  23^ 
and  BR+1  =  (2  +  l)BR  modulus  2& 
In  essence  they  have  used  one  generator  to  select  numbers  from  the 
other.  This  generator  seems  to  provide  an  improvement  in  some 
local  randomness  properties.  However,  the  penalty  for  the  improve- 
ment is  doubling  the  time  of  generation.  Marsaglia  and  MacLaren 
also  noted  that  the  method  of  table  look-up  may  once  more  become 
feasible.  In  the  case  of  the  CDC  I6OI4.  this  method  is  not  practical, 
However,  in  a  parallel  program  c  omputer  a  method  using  a  pair  of 
generators  may  well  be  advantageous.  The  generators  continually 
fill  up  the  bottom  of  a  short  table  in  memory  as  the  main  program 
uses  numbers  from  the  top.  The  size  of  the  table  is  chosen  to 
ensure  that  the  program  never  uses  all  the  numbers  in  the  table 
and  thus  the  effective  generation  time  will  be  just  the  load  cycle 
time. 

3.2.6  The  CODAP,  Fortran  63,  machine  language  program  for  the 
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Cumulative  Mean  Of  the  First  I1  Samples 
Of  The  Uniform  Random  Number  Generator 
(with  1000  numbers  per  sample) 
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Scatter  Diagram 
For  The  Uniform  (0,1)  Random  Number  Generator 
(Scale:  0.2  units  per  inch) 
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uniform  generator  used  as  a  basis  for  this  thesis  is  in  Appendix 
I.  The  expected  time  of  generation  per  number,  as  calculated  over 
several  samples  of  varying  size  was  found  to  be  552  microseconds 
per  number.  This  amounts  to  producing  1811  numbers  per  second. 
However,  the  generator  is  theoretically  much  faster  than  that.  The 
time  per  number  as  calculated  from  times  in  the  Control  Data 
Corporation  specifications  for  the  computer  is  121  microseconds. 
Measured  times,  depending  on  the  context  and  the  timing  mechanism 
varied  from  a  minimum  of  370  to  a  maximum  of  700  microseconds. 
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U.     The  Normal  Disbribution  (Univariate  Case). 
U.l    Distribution  characteristics. 

It  is  desired  to  generate  numbers  such  that  the 
density  function  will  be 

where  u  is  the  mean  of  the  distribution  and  o2  is  the  variance.  For 

the  basic  case  the  mean  is  taken  to  be  zero  and  the  variance  to  be 

one.  The  first  four  moments  are  Ml  =  0,  M2  -1,  M3-0,  MU-3. 

If  the  desired  distribution  is  to  have  a  mean  other  than  zero, 

say  u,  and  a  variance  other  than  one,  say  V,  then  the  following 

transformation  applied  to  the  numbers  generated  by  the  N(0,1) 

generator  developed  here,  represented  by  RN01,  will  produce  a 

number,  RNUV,  With  the  desired  characteristics. 

RNUV  =  (RN01)  (V)  +  U 

In  this  paper  the  normal  distribution  is  treated  in  three 

separate  sections.  The  univariate  case  is  developed  first,  then 

the  bivariate,  and  finally  a  general  multivariate  case  is  demonstrated. 

The  main  purpose  of  this  separate  treatment  is  to  allow  a  more 

efficient  handling  of  the  more  commonly  used  univariate  and 

bivariate  cases.  A  general  n-dimensional  normal  random  number 

generator  would  be  much  slower,  when  used  for  n  equals  one,  than 

the  univariate  generator  demonstrated  in  Section  U.  The  test 

procedures  for  each  generator  are  also  different. 

U.2    Methods  of  generation. 

U.2.1  The  normal  distribution  is  one  of  the  most  used  and 

15 


tabulated  distributions.  As  with  the  uniform  distribution  several 
procedures  have  been  used  to  produce  random  numbers  from  it.  The 
methods  discussed  here  are  those  that  are  most  adaptable  to  use 
on  a  computer* 

U.2.2  The  most  common  procedure  has  been  to  take  the  sum  of  K 
uniform  random  numbers.  The  central  limit  theorem  shows  that  this 
(with  the  mean  subtracted,  and  divided  by  the  standard  deviation) 
approaches  the  normal  as  K  gets  large.  Vaa  tested  a  generator 
using  the  sum  of  twelve  uniform  random  numbers.  This  approximation 
has  the  disadvantage  of  being  truncated  at  plus  and  minus  six. 
Even  more  important  a  factor  is  the  time  required  to  generate  these 
numbers .  It  is  hoped  that  a  more  exact  and  faster  method  can  be 
found. 

U.2.2  The  so-called  half -Gaussian  method  [ll]  provides  a  theoretically- 
exact  transformation  from  the  uniform  to  the  normal.  However  in  an 
attempt  to  reduce  time  some  fairly  drastic  approximations  have  been 
made.  These  approximations  should  not  affect  randomness,  however, 
and  should  only  be  a  factor  in  accuracy  beyond  the  fifth  decimal 
place.  The  following  flow  chart  is  the  basis  for  the  routine.  R(J) 
is  a  uniform  random  number.  HGRN  is  the  normal  number » 
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GENERATE  R(J),  R(J+l) 


Y 

Z 


-LN(R(J)) 
-LN(R(J+1)) 


J  =  J+2 


The  routine  first  generates  the  positive  half  of  a  normal  distribution 

then  adds  a  random  sign  selected  by  another  uniform  number*  This 

program  makes  no  external  calls  to  the  log  function  but  rather 

uses  the  following  series  approximation! 

T  =   (R-1)/(R+1) 
LN(R)=  2(T+t3/3+  TV5+  T7/7+  T9/9) 

The  CODAP  function  sub-program  is  Appendix  II. 

U.2.3  An  excellent  approximation  technique  has  been  developed  by 

Marsaglia  and  Bray  (J?J  .  Marsaglia  has  developed  several  similar 

techniques  but  the  one  proposed  here  seems  optimum  in  terms  of  time 

required  per  number  and  the  storage  space  required.  The  method 

involves  selecting  one  of  four  functions  of  varying  complexities  to 

produce  the  random  number.  One  function  is  very  simple  and  fast  and 
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is  used  86$  of  the  time;  the  next  function  is  also  fast  and  is 
used  11$  of  the  time.  The  remaining  three  percent  of  the  time 
much  more  complex  functions  are  used,  however,  due  to  the  rapidity 
with  which  97$  of  the  numbers  are  formed,  the  overall  expected 
generation  time  per  number  is  relatively  low.  The  program  outline 
is  as  follows  (MSRN  is  the  desired  random  number): 
1.  86.38$  of  the  time,  set 

MSRN  =  2[R(j)  +  R(J+1)  +  R(J+2)  -  1.5] 
2»  11.07$  of  the  tijne,  set 

MSRN  =  l.,5(R(J)  ♦  R(J+1)  -  l| 
3.  2.28002039$  of  the  time  form  pairs  (X,Y)  such  that 
X  =  6R(J)  -  3   and 
Y  =  0„358  R(J+1) 

until  Y*G3(X);  then  set  MSRN  =  X,  G3(X)  is  defined  by: 
17.U9731196exp(-xV2)-U.73570326(3-xi)-2a57875Uli(1.5-lx|)  for  bcKl 
17.u9731196exp(-xV2)-2.36785l63(3-lxlf  -2.157875UU(1.5-Ix|)  for  Kx<1.5 
17.u9731196exp(-xV2)-2.3678£L63(3-lxO  for  1.5<x<3 

U.  0.26997961$  of  the  time  form  pairs  (X,Y)  until  either  X 
or  Y  is  greater  than  three,  then  let  that  one  equal  MSRN.  For 
RM(J)  uniform  on  the  interval  (-1,1)  and  such  that  if  RM(J)X+  RM(J+l)*£l, 
then  let  Z  =  RM(J)*+  RM(J+lf 
and  X  =  RM(J)[  f 9-2LN(Z)j/(Z)] 

Y  =  RM(J+1)[  (9-2LN(Z)J/(Z)] 
The  CODAP  function  sub -program  is  Appendix  III. 
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U.3    Selection  of  generators, 

A  complete  table  of  results  of  the  tests  on  the  normal 
generators  is  in  Appendix  IV.  Partial  results  are  presented  below. 


Average  observed 
generation  time 
per  number 

Sample  range 

Sample  size 

Mean 
(theor.  =  0) 

Variance 
(theor.  *  1) 

3rd  Mom. 
(theor.  =  0) 

Uth  Mom. 
(theor,  =3) 

JN  X 
(K#05J..6U) 

fUftf-1 1 
J20WL) 

DmaxyN 
1.22//N 


HALF -GASS IAN  TECHNIQUE 
3625  microseconds 


MARSAGLIA  TECHNIQUE 
1503  microseconds 


1-1000       1-10000       10001-20000       1-1000       1-10000       10001- 

20000 


10000  10000  10000 

0.00  -0.02  -0.01 

.99  1.00  .99 

.01  -0.0U  -0.07 

2.77  2.97  2.93 

0.13  -2.U7  -1.06 

-O.lii  -0.11  -0.80 

0.0082  0.0095  0.00U9 

0.0386  0.0122  0.0122 


1000 

10000 

10000 

o.ou 

-0.01 

0.00 

1.06 

1.0U 

1.01 

-0.05 

-0.02 

0.01 

3.17 

3.20 

3.12 

1.32 

-0.96 

0.30 

1.30  2.57  0.95 
0.0270  0,0073  0.0033 
0.0386    0.0122    0.0122 


The  samples  generated  by  the  half -Gaussian  technique  passed  the  hypothesis 
on  the  mean,  at  the  10$  level,  80$  of  the  time ;   on  the  variance  at 
the  10$  level,  30$  of  the  time.  At  the  5$  level  the  test  on  the 
mean  was  passed  90$  of  the  time,  the  test  on  the  variance  was  passed 
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$0$  of  the  time.  The  Marsaglia  technique  generated  samples 
that,  at  the  10$  level,  passed  the  test  on  the  mean  80$  of 
the  time,  and  on  the  variance  90$  of  the  time.  The  Marsaglia 
technique  consistently  passed  the  Kolmogorov-Smirnov  test  but 
appeared  a  little  heavy  in  the  'tails1  none  the  less.  For  a 
further  examination  of  this  see  Addendum  1.  The  half -Gaussian 
technique  failed  the  Kolmogorov-Smirnov  test  once  but  appeared 
better  behaved  in  the  tails.  The  graph  on  page  21  is  a  scatter 
diagram  for  the  numbers  produced  by  the  Marsaglia  technique. 
The  graph  consists  of  U£00  points.  Marsaglia' s  technique  produced 
better  results  for  the  first  four  moments.  The  decisive  factor 
that  leads  to  the  selection  of  one  of  the  generators  is  the  time 
to  generate  each  number.  The  Marsaglia  technique  is  more  than 
twice  as  fast  as  the  half -Gaussian  technique,  is  the  one  used  in 
further  developments,  and  is  the  one  recommended  for  general  use» 
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Scatter  Diagram 
For  The  Normal  Random  Number  Qenerator 
(Scales  1.0  units  per  inch) 
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5.     The  Normal  Disbribution  (Bivariate  Case) 
5.1    Disbribution  characteristics. 

It  is  desired  to  generate  pairs  of  numbers  (X,(  >XZ"L  ) 
such  that  the  two-dimensional  random  variable  (l^t   X2^  has  tlie 
joint  density  function 

for  all  (x-ijXp).  The  constants  are  u^  and  iu,  the  means ; Cr%  (>  0),  <TX 
(>()),  the  standard  deviations;  and  £>  (-l^££l),  the  correlation 
coefficient.  Thus  the  requirement  for  a  random  vector  must  be 
accompanied  by  the  specification  of  the  mean  Vector  (u^Up),  the 
variances  (the  squares  of  the  standard  deviation),  and  the 
correlation  coefficient.  Another  equivalent  form  of  the  input 
would  be  the  mean  vector  and  the  covariance  matrix.  This  last 
form  of  input  will  be  used  in  the  general  multivariate  case. 
The  distribution  is  specified  in  matrix  notation  as  follows: 

f(i)  -  fOd,^)  -  lfi£«p£-±(Y-U>'R(r-u)J 

(2ry 

where  f(X)  is  the  joint  density  function  of  the  x'sj  p  is  the 
dimension-in  this  case  two;  U  is  the  mean  vector^  and  R  is  the 
inverse  of  the  convariance  matrix.  Thus  IRI   is  the  square  root 
of  the  determinant  of  the  inverse  of  the  covariance  matrix;  (Y-U) 
'R(Y-U)  is  a  quadratic  form,  where  (Y-U)1  is  a  one  by  p  vector, 
R  is  a  p  by  p  matrix,  and  (Y-U)  is  a  p  by  one  vector. 
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5.2    Method  of  generation. 

If  R(J)  is  a  normal  random  number,  the  random  bivariate 
vector  (VN1,  VN2)  is  formed  as  follows: 

VN1  -  (cr,  )R(J)+u, 


VN2  =  (ot)R(J)+  (  Pi)R(J+l)(JT^)+U2 
The  source  of  the  normal  random  numbers  is  the  Marsaglia  routine 
described  in  the  previous  section.  The  generator  is  Appendix  V. 
5.3    Testing  the  generator. 

First  the  maximum  likelihood  estimators  of  the  mean  vector 
and  the  covariance  matrix  are  formed  fl2]  .  The  maximum  likelihood 
estimators  for  the  parameters  ares 

Sit  =  JJP  I$?  (§,2  . 


The  distribution  of  the  mean  when  the  covariance  matrix  is  unknown 
was  shown  by  Hotelling  to  be  a  multivariate  analogue  of  the  t-test 
and  is  called  the  generalized  T  statistic.  However,  the  covariance 
matrix  is  known  and  once  again  we  can  use  a  more  powerful  test. 

H0  :  u  =  (U10  U20)«        HAs  u*  (Ui0  U20) ' 
Construct  the  statistic  H  such  that: 

H  =  N(xi-u10,  x2-u20)C"1(x1-u10,  x2-u20)' 
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where  C  is  the  given  covariance  matrix.  If  H  <  X2  0*)we  accept  the 
null  hypothesis.  The  test  for  a  hypothesized  mean  vector 
U  =  (0  ,  0)',  and  a  covariance  matrix  with  variances  one  and 
correlation  coefficient^,  reduces  to  calculating  H  such  that: 
N   /  7     T\    I  zr,   "  =§*  \  (  Xx 


-  ±L    (Xf-Zett   +-X2*  ),  for     cr.rcr-.r/. 

for  sample  sizes  of  1000  the  test  will  be  made  for^  =  0.25,  0.5, 

0.75.  At  the  10$  level  if  H  —  U.6l  accept  the  hypothesis.  In 

order  to  test  whether  the  covariance  matrix  is  a  given  matrix, 

the  test  as  outlined  in  Anderson  [12]   could  be  used,  but  before 

this  a  digression  into  the  philosophy  of  testing  is  in  order. 

The  purpose  of  these  hypothesis  tests  is  in  general  to  render  a 

judgement  as  to  whether  a  sample  of  data  from  some  type  of 

experiment  is  distributed  according  to  a  certain  distribution 

with  certain  parameters.  This  type  of  t est  was  applicable  to  our 

uniform  random  number  generator,  but  its  use  seems  extraneous  for 

the  distributions  derived  from  this  generator.  The  trans - 

form  itions  from  the  uniform  to  the  normal  to  bivariate  normal  are 

either  theoretically  exact  or  controllably  close  to  being  exact. 

What  is  really  needed  is  a  method  to  ensure  that  any  'inaccuracies 

there  may  have  been  in  the  basic  generator  are  not  somehow 

amplified  in  the  transformation  so  as  to  bias  the  derived  distribution. 
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To  this  end  sophisticated  testing  procedures  do  not  seem  to  be 
in  order.  Instead  the  testing  will  be  restricted  to  calculation 
of  estimators,  and  some  goodness  of  fit  tests.  The  testing  plan 
for  the  bivariate  generator  is  to  test  10  sets  of  1000  each  for 
each  of  several  correlation  coefficients.  The  generator  was  also 
tested  for  some  odd  combinations  of  parameters-  such  as  U=(100.0 
-100.0),  ^=0.$,  or^25.0,(g=  -0.8. 

With  0.75  as  the  correlation  coefficient,  the  fourth  set  showed  a 
disagreeably  large  difference  in  the  maximum  likelihood  estimators, 
The  H  statistic  was  close  to  the  critical  value  at  the  10$  level. 
The  sample  passed  the  test  in  90$  of  the  cases.  With  0.50  as  the 
correlation  coefficient  the  sample  passed  the  H  test  90$  of  the 
time.  Sample  set  four  once  again  had  rather  low  covariance 
estimators  and  once  again  had  3.76  as  the  value  for  H  (compared 
with  U.6l  for  the  critical  value).  Sample  set  six  was  again  the 
culprit  in  not  passing  the  H  test.  A  similar  result  was  noted 
in  the  set  using  0.25  as  a  correlation  coefficient.  This  suggests 
the  advisability  of  further  testing  of  the  Marsaglia  generator. in 
these  ranges.  Complete  test  results  are  in  Appendix  VI.  The 
generator  performs  as  expected  and  is  recommended  for  use. 


25 


— 

Pll  0*  •    0    " 
P21  P22--0 

«                *              *  • 

"Pll  P21  '—ML 

0       P22  "■•PN2 

•• 

>                 •     * 

PN1  PN2     PMN 

i 

0  .•► o  PNN. 

6.     The  Normal  Disbribution  (Multivariate  Case)* 

6.1  Distribution  characteristics. 

As  is  noted  in  Section  6.1  the  desired  distribution  is 
f(X)  =  f(x()Xi,...,xH)  -  -j^l  &cf(-i(\-v)'R^-o)) 
where  R  is  the  inverse  of  the  covariance  matrix. 

6.2  Method  of  generation. 

As  shown  by  Wold  [l6]  ,  the  method  is  based  on  a  triangularization 
of  the  covariance  matrix  such  that  if: 


Cll  C12  .  i .  C1N 
C21  C22  -••  C2N 


CN1  CN2  "  CNrt 

then  the  N-dimensional  random  vector  Z  may  be  formed  as  follows. 
The  RN(l)  are  the  normal  random  numbers  generated  by  the  Marsaglia 
technique. 

Z(l)  =  PHx  RN(1) 
Z(2)  =  P21  x  RN(1)  +  P22  x  RN(2) 
Z(3)  -  P31  x  RN(1)  +  P32  x  RN(2)  +  P33  x  RN(3) 
Z(N)  =  PN1  x  RN(1)  +  PN2  x  RN(2)  ♦  ...  +  PNN  x  RN(N). 
The  method  as  programmed,  assumes  all  the  means  are  zero,  but  simple 
addition  of  the  mean  when  required  will  remedy  this.  The  matrix 
triangularization  is  based  on  a  symmetric,  positive  definite  or 
positive  semi -definite  matrix  C.  Thus  any  pair  of  the  random 
variables  can  have  a  partial  correlation  coefficient  of  one.  The 
routine  will  set  all  elements  of  the  vector  that  are  dependent  equal 
to  zero.  This  procedure  was  selected  since  in  order  to  relate  the 


26 


variables  properly  requires  an  inordinately  extra  amount  of 
computer  time-  time  that  even  when  not  needed  adds  to  execution 
time.  If  the  user  desires  some  variable  to  be  a  linear  trans- 
formation of  some  others,  then  he  only  needs  to  keep  track  of 
which  variables  these  are,  and  where  the  program  sets  the  vector 
element  equal  to  zero,  substitute  the  appropriate  linear 
combination  of  the  corresponding  independent  elements  of  the 
vector.  The  routine  also  checks  and  where  the  dependence  is  very 
close  to  one  will  assume  a  correlation  of  one,  and  proceed  as 
noted  above.  This  procedure  is  required  to  prevent  division  by 
numbers  very  close  to  zero.  The  following  example  will  clarify 
the  above  explanation.  Suppose  it  is  desired  to  generate  random 
vectors  from  a  distribution  with  mean  vector  one  and  covariance 
matrix  C,  where: 

2  1 
C  -  1  It 

3  5 

Thus  column  three  is  a  linear  combination,  in  fact  the  sum.of 
columns  one  and  two.  This  means  that  the  third  variable  is  not 

an  independent  variable.  The  triangulation  routine  will  produce  a 
matrix  P  of  the  form: 


2     1 

3 

1     k 

5 

3    5 

8 

p  = 


a  o  o 
b  d  o 


£  e  o 
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As  expected  column  three  is  identically  zero.  Thus  if  the  first 
two  normal  (0,1)  numbers  generated  are  denoted  by  Rl  and  R2,  the 
generated  vector  will  be  of  the  form: 

Z  =  (aRl,   bRl  +  dR2,   0)   . 
Since  it  has  been  determined  that  the  third  variable  is  the  sum 
of  the  first  two  and  also  that  we  desire  all  the  means  to  be  one 
then  the  desired  vector  is  of  the  form: 

Z  =  (aRl  +  1,   bRl  +  dR2  +  1,   aRl  +  bRl  +  dR2  +  1) 
The  generator  is  Appendix  VII. 
6.3    Testing  the  generator. 

The  maximum  likelihood  estimators  of  the  mean  vector  and 
the  covariance  matrix  are  formed  as  follows: 

BMl(j)  =  jrQ   Y(I.«     r°r  T-l^...,* 

C(I,J)  =  ^[|  (Y(I,KW(T,«)}]   forr,T=«,2,.^. 

N  is  the  dimension  of  the  covariance  matrix,  M  is  the  number  of 

sample  vectors,  the  l(I,J)  are  the  vector  elements,  the  BMl(j) 

are  the  elements  of  the  mean  vector  estimates,  and  the  C(I,J)  are 

the  elements  of  the  covariance  matrix  estimate.  The  results  of 

the  tests  for  various  covariance  matrices  are  listed  in  Appendix  VIII, 
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7.  The  Poisson  Distribution. 

7.1         Distribution  characteristics. 

It  is  desired  to  generate  numbers  such  that: 
P  [x  *  a]    ■     e  "*m  ma       for  m>o;aasal,    ... 


(X 


Thus  the  Poisson  distribution  is  a  discrete  distribution  for 

integer  values  of  a,  and  is  characterized  by  the  parameter  m. 

2 
The  first  four  moments  are  Ml  =  M2  ■  M3  ■  m,  and  MiU  =  3m  +m. 

7.2    Method  of  generation. 

This  method  is  due  to  Kahn  [111  an<^  is  a  theoretically- 
exact  transformation.  The  flow  chart  of  the  routine  is  shown  belowj 


K  »  2 
Y(l)   =  1 

1 

t 

^ 

Y(K)   =  R(K-l)xY(K-l) 

\ 

i 

NO 

^^    lb 

K  ■  K+i 

<y(k)£  r*  ?  ^ 

YES 


POIS  is  the  generated  Poisson  random  number,  m  is  the  parameter,  E 
is  the  irrational  number  2#7l828.«^ 
7 #3    Testing  the  generator* 

The  first  four  moments  are  calculated  for  samples  of  $000 
or  1000  with  the  parameter  m  taking  on  various  values.  The 
Kolmogorov-Smirnov  goodness  of  fit  test  is  applied  to  some  of  the 
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samples.  As  Tate  and  Clelland  [l9j  have  stated  the  test  is 
applicable  to  discrete  distributions  with  neglible  changes  in 
significance  for  large  sample  sizes.  The  generator  performed 
well  and  consistently  passed  the  test.  The  test  results  are 
in  Appendix  X,  while  the  generator  itself  is  in  Appendix  IX. 
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8.     The  Exponential  Disbribution. 

8.1  Distribution  characteristics . 

It  is  desired  to  generate  random  numbers  such  that  the 
density  function  is:  f(a)  =  e"a. 

The  first  four  moments  are  Ml  =  M2  =  1,  M3  =  2,  Ml*  ■  9-  The 
distribution  is  often  specified  as: 

f(a)  =  X  e"**"    for  A>oya*o 

where  A  is  the  parameter  of  the  distribution.  The  generator  here 
takes  the  case  where  A =  1.  However  the  exponential  distribution 

has  the  characteristic  that  if  the  numbers  generated  here,  (Expl) 
for  which  A  =  1,  are  simply  multiplied  by  the  parameter  desired 
for  the  distribution  (CAMBDA),  then  the  desired  numbers  are  generated. 
EXPL  =  EXP1*  LAMBDA 

8.2  Method  of  generation. 

This  method  is  a  theoretically  exact  procedure  due  to 
Marsaglia  [13]  •  A  more  obvious  method  would  be  to  take  the 
integral  transformation-  the  negative  logarithm  of  a  uniform 
random  number.  However  this  method  is  slower  on  most  computers 
than  the  one  demonstrated  here.   Let  C  =  l/(e-l),  and  let  the 
random  variable  N  take  on  the  values  1,2,3,U,...  with  probabilities 
c,  c/2'.  ,0/3  I  ,...  Then  let  the  random  variable  M  take  values  0,1, 
2,3,...  with  probabilities  1/ce,  1/ce  ,  1/ce^,...  Then  we  form  the 
desired  random  number: 

EXP1  =  M+MIN(U1,  U2,...,UN) 

The  CODAP  generating  routine  is  Appendix  XI.  The  sample  moments  and 
the  results  of  the  goodness  of  fit  tests  are  in  Appendix  XII. 

31 


9.     A  General  Disbribution. 

The  method  of  obtaining  random  numbers  described  in  this 
section  is  applicable  only  to  a  restricted  set  of  distributions. 
However,  the  method  is  applicable  to  a  type  of  distribution  that 
frequently  occurs  in  model  building  and  war  gaming.  The  method 
is  examined  by  use  of  an  example.  A  discussion  of  how  and  when 
to  use  this  method  for  other  distributions  is  included. 

It  is  desired  to  produce  numbers  from  the  density  function 
f(x)  such  that: 

f(x)  =  2-2x   for  Otx£l 
The  method  is  based  on  drawing  uniform  numbers  in  pairs,  normalizing 
the  scale  of  the  numbers,  and  testing  to  see  if  the  point  formed 
by  the  pair  lies  under  the  curve  described  by  the  density  function. 
If  it  does,  the  x  coordinate  of  the  point  is  taken  as  the  random 
number;  if  it  does  not,  another  pair  of  uniform  numbers  are  drawn 
and  the  process  is  continued.  The  flow  chart  for  f(x)  =  2-2x  is: 


GENERATE 
U(J),  U(J+1) 

t 

t 

Y  =  2(U(J+I)) 

i 

^^      IS    ^^^ 

NO 

T     a      Tj.O 

^^^Y  =  2-2X  7^* 

u      m 

YES 


The  U(J)  are  the  uniform  (0,1)  random  numbers  and  RN  is  the  random 
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number  from  the  distribution  f(X).  Thus  the  method  could  be 
applied  with  theoretical  exactness  to  any  bounded  continuous 
distribution.  Another  commonly  used  density  is  g(x)  where 

g(x)  ■  nx        for  0<x<l,  nil 
This  density  function  is  bounded  and  continuous  and  a  member  of 
the  class  to  which  this  method  is  applicable.  In  many  applications 
the  user  must  make  a  judgement  about  the  amount  of  use  a  generator 
is  going  to  get.  If  it  is  intended  for  heavy  use  it  may  be  well 
to  explore  the  literature  for,  or  to  design,  a  more  efficient 
generator  than  the  type  described  in  this  section.  Hoxrever,  if  the 
generator  is  to  only  be  used  a  limited  amount,  or  if  only  a 
restricted  amount  of  resources  are  available,  this  generator  is 
easily  programmed  and  will  be  eminently  satisfactory  in  a  wide 
variety  of  cases.  The  efficiency  of  the  generator  may  be  defined 
as  the  reciprocal  of  the  expected  number  of  iterations  needed  to 
produce  a  point  under  the  curve. 


m 


^-*x 


>K 


»  X 


la.  ">.  le. 

Figure  la  clearly  has  an  efficiency  of  1/2,  figure  lb  has  an 

efficiency  of  less  than  1/2,  and  as  n  gets  large  the  efficiency 
decreases  rapidly. 
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A  distribution  such  as  figure  lc  where  h(x)  goes  to  zero  at  some 
large  value  of  x  may  vary  well  have  such  a  low  efficiency  as  to 
make  this  generator  unsuitable.  If  h(x)  only  approaches  zero  as 
x  becomes  large,  and  therefore  x  is  not  bounded,  then  another 
generator  should  be  used.  However,  if  one  cannot  be  found  and 
the  user  is  not  too  concerned  with  the  tail  of  the  distribution, 
we  can  easily  adapt  this  generator.  Suppose  we  wished  to  generate 
numbers  for  the  Weibull  distribution,  which  is  a  three  parameter 
distribution  used  widely  in  reliability  theory,  then  the  following 
procedure  might  be  followed: 


Since  the  distribution  may  take  any  of  the  forms  in  figure  2,  we 
must  first  limit  the  development  to  those  parameter  combinations 
that  are  bounded  at  x  =  0.  Since  the  user  is  also  often  concerned 
with  the  behaviour  in  part  of  the  tail  it  must  first  be  decided  if 
the  distribution  could  simply  be  truncated  at  some  point.  Even  if 
this  is  acceptable  an  efficiency  of  .00001  can  easily  be  envisioned. 
If  it  is  unacceptable,  the  tail  beyond  some  point  could  be  closely 
approximated  by  the  exponential  distribution. 
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10.    Summary  and  Conclusions. 

The  analyst  who  desires  to  use  any  of  the  generators 
demonstrated  here  is  once  again  warned  that  although  the  generators 
are  recommended  for  general  use  they  all  have  aberrations  of  some 
kind.  The  general  form  of  these  aberrations  is  noted  where 
known.  If  a  user  suspects  from  analysis  of  his  results  that  the 
aberration  in  the  generator  is  influencing  his  results  then  it  is 
recommended  that  he  first  modify  the  uniform  generator  by  one  of 
the  methods  suggested  in  section  two.  The  chance  of  this  occuring 
is  remote  and  other  parts  of  the  model  should  be  checked  carefully. 

The  programs  demonstrated  here  are  very  fast.  The  problem 
of  speed  is  stressed  here  and  is  a  major  decision  criterion  in  the 
selection  of  generators.  In  some  applications  the  speed  may  not 
be  as  important  a  factor,  and  in  this  case  the  type  of  generator 
discussed  in  section  nine  may  be  very  cost  effective  in  that 
little  investment  in  programming  and  testing  is  required.  The 
problem  of  measuring  the  speed  of  generation  of  a  number  is  not  as 
simple  as  it  may  appear.  The  generation  time  depends  on  the  program 
using  the  generator  and  also  on  the  method  used  to  time  the  routine. 
The  time  to  generate  a  number  based  on  the  Control  Data  Corporation 
specifications  for  the  I6OI4.  computer  theoretically  should  be  121 
microseconds.  The  observed  generation  times  vary  from  370  to  700 
microseconds. 

The  tests  used  here  are  statistically  sound,  but  the  meaning 

of  the  results  of  several  tests  of  the  same  generated  sample  could 
bear  further  investigation* 
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Some  very  interesting  new  methods  for  generating  random 
numbers  from  various  distributions  have  been  developed  by 
H.  Rubin.  Rubin's  work  is  soon  to  appear  as  a  Stanford  Applied 
Statistics  Laboratory  Technical  Report..  It  would  be  of  interest 
to  compare  his  methods  with  the  technique  used  above. 
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Appendix  I. 

The  Fortran  63  CODAP  function  subprogram  for  generating  uniform 

(0,1)  random  numbers. 

The  program  is  called  by  using  a  variable,  'UNIFORM(DUMMY) ' „  The 

argument  'DUMMY1  is  not  used.  For  example:  A  =  UNIFORM ( DUMMY )+  B. 

No  common  or  dimension  entries  are  required* 

The  observed  average  length  of  the  program  is  552  microseconds. 

The  starting  number  may  be  change  to  'R   OCT    5U02033U50727U22 ' 

if  it  is  desired  to  enter  the  generator  near  the  middle  of  the 

period. 


IDENT 

UNIFORM 

ENTRY 

UNIFORM 

Ml 

OCT 

U000000000000000 

M2 

OCT 

2000000000000000 

R 

OCT 

1777777777777777 

UNIFORM 

SLJ 

■H-* 

LDA 

■* 

ALS 

21 

INA 

1 

SAU 

EXIT 

ENQ 

0 

LDA 

R 

LLS 

+19 

SCL 

Ml 

ADD 

R 

ADD 

R 

ADD 

R 

STA 

R 

ARS 

11 

ADD 

M2 

FAD 

M2 

EXIT 

SLJ 

•*# 

-  Cf*- ) 


END 


39 


Appendix  II. 

The  half -Gaussian  technique  for  producing  normal  random  numbers. 

The  number  is  generated  by  the  use  of  the  expression  'GSRN(BLK)1 , 

The  number  is  returned  as  GSRN..  The  argument  .'BLK'  is  not  used. 

No  external  calls  are  made. 

The  observed  average  length  of  time  to  generate  one  number  is 

3625  microseconds. 


IDENT 

GSRN 

ENTRY 

GSRN 

GSRN    SLJ 

■M-* 

LDA 

♦ 

ALS 

2U 

INA 

1 

SAU 

GSRN 

.190    RTJ 

UNIFORM 

STA 

Bl 

FAD 

=02  001U00000000000 

STA 

TEMP 

LDA 

Bl 

FSB 

=02  001U00000000000 

FDV 

TEMP 

STA 

X 

FMU 

X 

STA 

X2 

FMU 

=0177U707070705726 

FAD 

=0l775UUiUiiUUi3O23 

FMU 

X2 

FAD 

-01775631U631U6315 

FMU 

X2 

FAD 

=017765252 52 52U3U1 

FMU 

X2 

FAD 

=02  001U 00000000000 

FMU 

X 

FMU 

=05775377777777777 

STA 

Y 

RTJ 

UNIFORM 

STA 

Bl 

FAD 

=02  OOLU  00000000000 

STA 

TEMP 

LDA 

Bl 

FSB 

=02  001U00000000000 

FDV 

TEMP 

Uo 


STA 

X 

FMU 

X 

STA 

X2 

FMU 

-0177U707  070705726 

FAD 

=01775UUiiUUUUi3023 

FMU 

X2 

FAD 

=0177563H|631U6315 

FMU 

X2 

FAD 

-01776525252 52U3U1 

FMU 

X2 

FAD 

=02  ooiu 00000000000 

FMU 

X 

FMU 

=0577U377777777777 

STA 

Z 

LDA 

Y 

FSB 

=02  ooiuooooooooooo 

STA 

YM1 

FMU 

no. 

THS 

z 

SLJ 

.190 

RTJ 

UNIFORM 

TK3 

=02 ooou 00000000000 

SLJ 

#+2 

LDA 

Y 

SLJ 

GSRN 

LAC 

Y 

UNIFORM 

SLJ 

GSRN 

SLJ 

#* 

ENQ 

0 

LDA 

R 

LLS 

19 

SLC 

=0U 000000000000000 

ADD 

R 

ADD 

R 

ADD 

R 

STA 

R 

ARS 

11 

ADD 

=02000000000000000 

FAD 

=02000000000000000 

SLJ 

UNIFORM 

R 

OCT 

177777777777777777 

Bl 

BSS 

1 

TEMP 

BSS 

1 

X 

BSS 

1 

Y 

BSS 

1 

YM1 

BSS 

1 

Z 

BSS 

1 

X2 

BSS 

1 

X3 

BSS 

1 

X5 

BSS 

1 

X7 

BSS 

1 

X9 

BSS 
END 

1 

Ul 


Appendix  III, 

The  number  is  generated  by  use  of  the  variable  •MARS(ZQ)1.  The 

argument  »ZQ'  is  not  used.  For  example:   'A  «  MARS(ZQ)/9. 

External  calls  are  made  to  LOGF,  SQRTF,  and  EXPF. 

No  common  dimension  entries  are  required. 

The  observed  average  length  of  time  to  generate  one  number  is 

1503  microseconds. 


IDENT 

ZMARS 

ENTRY 

ZMARS 

UNIFORM    SLJ 

** 

ENQ 

0 

IDA 

R 

LLS 

19 

SCL 

-ouooooooooooooooo 

ADD 

R 

ADD 

R 

ADD 

R 

STA 

R 

ARS 

11 

ADD 

=02  000000000000000 

FAD 

=02  000000000000000 

SLJ 

UNIFORM 

R         OCT 

1777777777777777 

ZMARS      SLJ 

*-* 

LDA 

* 

ALS 

2U 

INA 

1 

SAU 

ZMARS 

RTJ 

UNIFORM 

Gl        THS 

PI 

SLJ 

G2 

RTJ 

UNIFORM 

STA 

NORM 

RTJ 

UNIFORM 

FAD 

NORM 

STA 

NORM 

RTJ 

UNIFORM 

FAD 

NORM 

FSB 

=02001600000000000 

STA 

NORM 

FAD 

NORM 

SLJ 

ZMARS 

U2 


G2 


G3 

G31 


G3P0S 


T1.5 


THS 

P2 

SLJ 

G3 

RTJ 

UNIFORM 

STA 

NORM 

RTJ 

UNIFORM 

FAD 

NORM 

FSB 

-02  OOLUOOOOOOOOOOO 

FMU 

=02  001600000000000 

SLJ 

ZMARS 

THS 

P3 

SLJ 

GU 

RTJ 

UNIFORM 

FMU 

-02  00U600000000000 

FSB 

-02  002  600000000000 

STA 

T 

AJP 

2  GSFOS 

SCM 

•07777777777777777 

STA 

T+l 

FMU 

T+l 

STA 

T+2 

FMU 

=05777377777777777 

CALL 

EXPF 

FMU 

CI 

STA 

T+3 

LDA 

T+l 

THS 

=02 OOlU 00000000000 

SLJ 

T1.5 

LDA 

=02002600000000000 

FSB 

T+2 

FMU 

C2 

FAD 

T+3 

STA 

T+3 

LDA. 

=02001600000000000 

FSB 

T+l 

FMU 

C3 

FAD 

T+3 

SLJ 

G3END 

THS 

=02  001600000000000 

SLJ 

T3.0 

LDA 

=02002600000000000 

FSB 

T+l 

STA 

T+U 

FMU 

T+U 

FMU 

CU 

FAD 

T+3 

STA 

T+3 

LDA 

=02  001600000000000 

FSB 

T+l 

FMU 

C3 

FAD 

T+3 

SLJ 

G3END 

U3 


T3.0 


G3END 


GU 


LDA 

=02002600000000000 

FSB 

T+l 

STA 

T+U 

FMU 

T+U 

FMU 

CU 

FAD 

T+3 

STA 

T+U 

RTJ 

UNIFORM 

FMU 

CS 

THS 

T+U 

SLJ 

G31 

LDA 

T 

SLJ 

ZMARS 

RTJ 

UNIFORM 

FMU 

=02002UOOOOOOOOOOO 

FSB 

=02 OOlU ooooooooooo 

STA 

T 

FMU 

T 

STA 

T+l 

RTJ 

UNIFORM 

FMU 

=02002UOOOOOOOOOOO 

FSB 

=02  OOlU ooooooooooo 

STA 

T+2 

FMU 

T+2 

FAD 

T+l 

THS 

=02  OOlU ooooooooooo 

SLJ 

GU 

STA 

T+l 

CALL 

LOGF 

STA 

T+3 

FAD 

T+3 

SCM 

=07777777777777777 

FAD 

=02  OOUUUOOOOOOOOOO 

FDV 

T+l 

GALL 

SQRTF 

STA 

T+3 

FMU 

T 

STA 

T+U 

AJP 

2  *+l 

SCM 

=07777777777777777 

THS 

=02002600000000000 

SLJ 

GOOD 

LDA 

T+3 

FMU 

T+2 

STA 

T+U 

AJP 

2  *+l 

SCM 

=07777777777777777 

THS 

=02  002600000000000 

SLJ 

GOOD 

uu 


SLJ 

GU 

GOOD 

LDA 

T+U 

SLJ 

ZMAilS 

CI 

DEC 

17.U9731196 

C2 

DEC 

-U. 73570326 

C3 

DEC 

-2.157875UU 

CU 

DEC 

-2.36785163 

C5 

DEC 

0.358 

PI 

DEC 

O.8638 

P2 

DEC 

0.97U5 

P3 

DEC 

0.9973002039 

T 

BSS 

5 

NORM 

OCT 
END 

0 

U5 


Appendix  IV. 

Test  results  for  Marsaglia  normal  generator. 


The  test 

consists  of 

ten  consec 

sutive  sampl 

es  with  IOC 

100  numbers 

in  each. 

MEAN 

M2 

M3 

MU 

JNX 

S-(N-ll 

d/n 

.0 

1.0 

.0 

3.0 

-1.6U5 
to  1.6U5 

-1.6U5 
to  1.6U5 

-.02U7U 

.998UU 

-.0UU35 

2.96870 

-2.U737 

-.110U 

0.0095 

-.01062 

.98869 

-.07U87 

2.92989 

-1.0625 

-.7996 

0.00U9 

-.0056U 

.99388 

-.07973 

2.93538 

-  .56U0 

-.1030 

0.0027 

-.00873 

1.00U52 

-.09225 

2.99577 

-  .8729 

.3196 

0.0060 

-.00921 

.999U5 

-.OU77U 

3.03286 

-  .9209 

-.0392 

o.oouo 

-.00683 

.99709 

-.10706 

3.08U87 

-  .683U 

-.2057 

O.OOI48 

.00028 

.99963 

-.OU737 

3.03271 

.0277 

-.026U 

0.0030 

.00732 

1.01508 

-.OU782 

3.16U72 

.7317 

1.066U 

O.OOUl 

-.00U10 

1.02U5U 

-•05629 

3.09820 

-  .UioU 

1.7352 

0.0050 

-•01732 

.97325 

-.05067 

2.79351 

-1.7319 

-1.8911 

0.0050 
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Appendix  IV.  (Cont'd) 

Test  results  for  half -Gaussian  method. 


The  test 

consists  oi 

ten  consec 

:utive  sampli 

as  with  10000  numbers 

in  each 

MEAN 
.0 

M2 
1.0 

M3 
.0 

Mil 
3.0 

-1.6U5  to 
1.6U5 

S-(N-l) 

-1.6U5  to 
1.6U5 

d/n 

-.00961 

1.03629 

^,01728 

3.1977U 

-.96lii 

2.5661 

0.0073 

.00301 

1.03U1 8 

.00973 

3.12078 

-.3008 

.9529 

0.0033 

-.01107 

1.02600 

-.01559 

3.12095 

-1.1066 

1.8386 

0.0056 

.008UU 

1.0U205 

-.02587 

3.23017 

.8UU3 

2.9735 

0.0072 

.00230 

1.01U10 

-.OOU58 

3.13768 

.2295 

.9969 

0.00U1 

-.OOU75 

1.02696 

-.00209 

3.276U2 

-J*7U7 

I.906I1 

0.0036 

-.022U6 

1.0UU39 

,0l6i|2 

3.29951 

-2.2U59 

3.138U 

0.01U9 

-.00372 

1.02968 

-.01805 

3.13UU2 

-.3716 

2.0986 

0.0065 

-.0068U 

1. 0U 955 

.OU399 

3.28683 

-.68I1I 

3.5036 

0.0076 

.01761 

1.01793 

-.00633 

3.11080 

1.761U 

1.2680 

0.008U 
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Appendix  V. 

The  bivariate  normal  random  number  generator  (using  Marsaglia's 

technique). 

The  number  is  generated  by  a  call    'CALL  AKHN(VN1,  VN2)».     The 

bivariate  vector  is  returned  as  the  arguments  VN1,  VN2. 

Subroutine  AKHN  has  the  following  arguments  in  common:  SIGLSQ, 

SIG2SQ,  RHO,  Ul,  U2.     The  SIGISQ  are  the  desired  variances,  RHO 

is  the  correlation  coefficient,   and  Ul  and  U2  are  the  desired 

means.     Thus  besides  reading  in  these  values  in  the  main  program 

they  must  be  communicated  to  the   subroutine  by  a  common  statement 

such  as:    •COMONAlG]£Q/SIGlSQ/SIG2SQ/SIG2SQ/RHO/RHO/Ul/u'^2/U2'< 

External  calls  are  made  to  LOGF,  SQRTF,  and  EXPF. 

The  observed  average  generation  time  for  a  pair  of  numbers  is 

6I48O  microseconds. 


SIGISQ 

SIG2SQ 

RHO 

Ul 

U2 

AKHN 


IDENT 
BLOCK 


AKHN 

1 


COMMON  SIGISQ (1) 

BLOCK  1 

COMMON  SIG2SQ(1) 

BLOCK  1 

COMMON  RH0(1) 

BLOCK  1 

COMMON  Ul(l) 

BLOCK  1 

COMMON  U2(l) 

ENTRY  AKHN 

SLJ  ** 

LDA  * 

ALS  2U 

SAU  *+2 

INA  1 

SAU  EXIT 

LDA  ** 

SAL  VN2 


U8 


ALS 

2U 

SAL 

VN1 

+ 

RTJ 

ZMARS 

+ 

STA 

Al 

IDA 

SIG1SQ 

+ 

CALL 

SQRTF 

+ 

FMU 

Al 

STA 

A2 

VN1 

FAD 

Ul 

STA 

VN1 

LDA 

A2 

FMU 

RHO 

FAD 

U2 

STA 

A2 

LAC 

RHO 

FMU 

RHO 

FAD 

-02 OOlU OOOOOOOOOOO 

FMU 

SIG2SQ 

+ 

CALL 

SQRTF 

+ 

STA 

Al 

+ 

RTJ 

ZMARS 

+ 

FMU 

Al 

VN2 

FAD 

A2 

STA 

VN2 

EXIT 

SIxJ 

■** 

ZMARS 

SLJ 

■** 

RTJ 

UNIFORM 

Gl 

THS 

PI 

SLJ 

G2 

RTJ 
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Appendix  VI. 

Test  results  for  bivariate  normal  generator. 

The  test  consists  of  ten  consecutive  samples  of  1000  numbers  each  for 

each  of  three  different  correlation  coefficients. 
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COVARI/INCES       CORE.  . 
COEFF. 
1.0     1.0               1.0 

Means 

1.0     1.0 

MAXIMUM  LIKE. 
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Appendix  VII 

The  FORTRAN  programs  to  produce  multivariate  normal  random  vectors. 
The  first  entry  must  be  'CALL  TRIANG!.  This  produces  the 
triangularized  matrix,  and  allows  repetitive  calls  to  'CALL  MULTN(Z)1 
The  argument  'Z'  is  the  starting  address  of  the  random  vector. 
Several  other  entries  are  necessary.  The  dimension  of  the 
covariance  matrix  (and  the  desired  vectors)  is  set  equal  to  'NR'. 
The  desired  covariance  matrix  is  stored  as  a  matrix  called  'C ' • 
A  common  statement  with  'NR'  and  'C,  with  'C1  ap propriately 
dimensioned,  is  included.   • Z '  is  also  dimensioned.  The 
following  sample  program  will  read  in  a  matrix  'C1,  which  is  a 
5  by  5  matrix  punched  column  by  column  on  cards.  The  random 
vectors  are  stored  in  a  5  by  100  array.  Subroutine  'MULTN' 
and  subroutine  'TRIANG'  follow.  Function  'ZMARS',  the  normal 
random  number  generator,  from  Appendix  III  must  be  added. 
Subroutine  TRIANG  makes  external  calls  to  SQRTF. 

PROGRAM  EXAMPLE 

COMMON/NR/NR/C /C ( 5 , 5 ) /P/P (5,5) 

DIMENSION  Z(5),  ARRAY(100,5) 

NR-5 

READ  101,  ((C(I,J),I-1,NR),J-1,NR) 
101  FORMAT  (5F8.U) 

CALL  TRIANG 

DO  111  J»l,100 

CALL  MULTN (Z) 

DO  111  K-1,5 
111  ARRAY(J,K)=Z(K) 

STOP 

END 

For  a  three  by  three  matrix  subroutine  TRIANG  takes  10600 
microseconds,  and  each  vector  is  produced  in  9520  microseconds. 
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It  is  noted  that  since  the  programs  are  written  in  FORTRAN, 

they  are  not  in  any  sense  optimal. 

SUBROUTINE  TRIANG 

COMMON/NR/NR/C/C(3,3)/P/P(3,3) 

NC-1 

dq:5  I  ■  1,  B 

DQ  5  J  -  1,  NR 

5  p(i,j)  ■  o.o 

P(1,1)=SQRTF(C(1,1)) 

IF(P(1,1).LE.. 0001)771, 9 
9   DO  15  K»2,NR 
1$  P(K,1)«C(K,1)/P(1A) 

18     DO  83  JH=2,NR 
Q=0.0 
NC»NC+1 
NCM-NC-1 
NCP-NC+1 
DO  50  J-1,NCM 

50  Q=Q+P(NC,J)*P(NC,J) 
X=C(JH,JH)-Q 

IF  (X.LE.. 000005)  777,505 
505  P(JH,JH)=SQRTF(X) 

51  DO  55  K=NCP,NR 
S-0.0 

DO  53  JQ-1.NCM 
53  S»S+P(K,JQ)*P(NC,JQ) 
55  P(K,NC)-(C(K,NC)-S)/P(NC,NC) 
83  CONTINUE 

GO  TO  666 
771  DO  773  L-1,NR 
773  P(L,l)-0.0 

GO  TO  18 
777  DO  780  WH,NR 
780  P(L,JH)«0.0 

GO  TO  83 
666  CONTINUE 

RETURN 

END 
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SUBROUTINE  MULTN(Z) 

C OMMON/NR/NR/P/P( 3 , 3 ) 

DIMENSION  Z(10),RZM(10) 

DO  220  J-l,  NR 

Z(J)-0.0 
220  RZM(J)"ZMARS(DUM) 

DO  230  L-1,NR 

DO  230  M-1,NR 
230  Z(L)-Z(L)+P(L,M)*RZM(M) 

RETURN 

END 
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Appendix  VIII 

The  tests  are  for  a  sample  size  of  100  vectors. 

The  first  matrix  tested  was  a  10  by  10  identity  matrix.  The 

maximum  likelihood  estimate  of  the  covariance  matrix  was: 
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Covariance 
Matrix  Input 

Triangularized 
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Estimate 
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Appendix  IX. 

The  Poisson  distributed  random  number  generator. 

The  numbers  are  generated  by  an  initial  use  of  the  variable 

'NPOISSET (mean) ' ,  which  initializes  the  generator  for  the  desired 

value  of  the  mean.  This  is  followed  by  calls  to  'NPOIS(DUM)'  to 

actually  produce  the  random  numbers.  The  argument  'DUM1  is 

not  used.  For  example  for  a  desired  mean  of  3.0: 

X(l)  =  NP0ISSET(3.0) 

• 

X(J)  «  NPOIS(DUM) 
External  calls  are  made  to  ,EXPF. 

The  observed  average  generation  time  per  number  was  found  to  be 
approximately  representable  by:  TIME  -  600  +  630(MEAN). 
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Appendix  X. 

Test  results  for  the  Poisson  generator. 

The  test  consists  of  ten  consecutive  samples  of  either  1000  or 

5000  numbers  each-for  various  parameter  values. 
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PARAMETER 

Ml 
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5.0 

5.0 

5.0 

80.0 

.039  (10$) 

10.0 

9.9760 

9.6751 

10.0922 

276.3385 

.013       1000 

6935 

9.8280 

9.5560 

9.8887 

280.6076 

.030 

9.9070 

10.U328 

12.0257 

322.7622 

.027 

10.15U0 

9.8UU1 

9.U111 

287.9662 

.019 

10.1090 

10.8179 

10.537U 

3U7.2039 

.021 

9.9710 

9.9922 

8.U5U8 

305.978U 

.012 

10.0520 

10.2676 

11.U099 

298.631U 

.029 

10.1500 

10.3258 

11.2U6U 

307.3533 

.020 

9.99UO 

9.U59U 

7.103U 

265.9739 

.017 

10.03U0 

9.6UU5 

8.7856 

283.3U59 

.011 

Theor. 

10.0 

10.0 

10.0 

310.0 

.39  (10£) 
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PARAMETER 

mi 

M2 

M3 

MU 

d/n 

Sample 
Size 

$.0 

U.93UO 

U.97U2 

5.3U30 

78.80U5 

.016 

5000 

5.0U62 

5.0975 

U.9371 

80.2735 

.01U 

5.0108 

U.93U5 

U.6862 

7U.9158 

.005 

U.9912 

5.0533 

5.1257 

80.938U 

.00U 

5.0038 

5.0U12 

5.7602 

87.0035 

.003 

5.0126 

5.0367 

5.339U 

83.8973 

.005 

5.03U6 

U.98Uli 

U.8206 

80.172U 

.010 

5.0016 

5.01U2 

5.1778 

77.76U8 

.006 

5.021U 

U.9U19 

U.7638 

75.7139 

.006 

U.9598 

5.011*0 

U.8610 

77.982U 

.009 

Theor. 

5.0 

5.0 

5.0 

80.0 

.017 

(10%) 

9.99U8 

IO.0720 

10.U326 

303.7092 

.006 

5000 

10.0U02 

9.9338 

9.U176 

292.8100 

.011 

10.017U 

10.1211 

12.32U6 

335.7205 

.00U 

10.0260 

9.8721 

9.3898 

296.6628 

.016 

10.0852 

9.85U7 

9.8097 

302.7969 

.017 

10.0176 

9.8136 

9.8867 

29k. 72  93 

.011 

9.991U 

IO.1918 

11.8066 

323.726U 

.006 

Theor. 

10.0 

10.0 

10.0 

310.0 

.017 

(10$) 

6U 


Appendix  XI. 

Exponential  random  number  generator. 

The  numbers  are  generated  by  a  call  to  'EXPRN(DUM) ' .  The  argument 

'DUM'  is  not  used.  For  example:  Y  -  EXPRN(DUM)  +  U.2 

No  external  calls  are  made.  The  observed  average  generation  time 

for  one  number  is  2270  microseconds.  A  conversion  for  numbers 

with  parameter  other  than  one  is  given  on  page  31  in  Section  8.1. 


IDENT 

EXPRN 

ENTRY 

EXPRN 

p 

DEC 

1.00000000 

DEC 

.9999998U 

DEC 

.9999982U 

DEC 

.99998381 

DEC 

.99986831 

DEC 

.9990600U 

DEC 

.99U21023 

DEC 

.96996120 

DEC 

.87296508 

DEC 

.58197672 

Q 

DEC 

1.00000000 

DEC 

.99999989 

DEC 

.99999970 

DEC 

.99999917 

DEC 

.9999977U 

DEC 

.99999386 

DEC 

.99998330 

DEC 

.99995U60 

DEC 

.99987659 

DEC 

.99966U5U 

DEC 

.99908812 

DEC 

.99752125 

DEC 

.99326205 

DEC 

.98l68ii36 

DEC 

.95021293 

DEC 

.86U66U71 

DEC 

.63212055 

TABLE 

BSS 

20 

MIN 

DEC 

0.0 

R 

OCT 

5UO2033U50727U22 

UNIFORM 

SLJ 

** 

ENQ 

0 

LDA 

R 

LLS 

19 
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SCL 

=ou 00000000000000 

ADD 

R 

ADD 

R 

ADD 

R 

STA 

R 

ARS 

11 

ADD 

-0200000000000000 

FAD 

-0200000000000000 

SLJ 

UNIFORM 

EXPRN 

SLJ 

■** 

LDA 

* 

ALS 

2k 

INA 

1 

SAL 

ih+1 

SIL 

1  Ah 

SIU 

2  AU+1 

ENI 

1 

+10 

RTJ 

UNIFORM 

+ 

THS 

1 

P 

SLJ 

Aii 

+ 

ENI 

2 

-1 

Al 

INI 

2 

+1 

RTJ 

UNIFORM 

+ 

STA 

2 

TABLE 

+ 

ISK 

1 

+9 

SLJ 

Al 

+ 

STA 

MIN 

A2 

UP 

2 

A31 

SLJ 

A3 

A31 

LDA 

2 

TABLE 

THS 

MIN 

SLJ 

A2 

STA 

MIN 

SLJ 

A2 

A3 

LDA 

MIN 

+ 

ENI 

1 

+17 

RTJ 

UNIFORM 

♦ 

THS 

1  C 

SLJ 

til 

ENA 

1 

0 

INA 

-16 

SCM 

-07777777777777777 

ADD 

=02  Ohk  000000000000 

FAD 

=02  014*  000000000000 

Aii 

FAD 

MIN 

ENI 

1 

** 

ENI 

2 

#* 

SLJ 

*# 

END 
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Appendix  XII. 

Test  results  for  the  exponential  generator. 

The  test  consists  of  consecutive  samples  of  1000  each  for  lambda 


equal  one 

1  • 

SAMPLE 

Ml 

M2 

M3 

Mh 

d/n         Generation 
Time 

1000-1 

1.0056 

.95UU 

1.5731 

6.0113 

.015 

-2 

1.0602 

1.1981 

2.93UO 

15.0U79 

.020 

-3 

1.0606 

1.2050 

2.9526 

15.5327 

.037 

-U 

1.0272 

1.1211 

2.5961 

12.1087 

.020 

-5 

.9122 

.8821 

1.7657 

7.6101 

.030 

-6 

1.02l4i 

.9U8U 

1.U382 

5.209U 

.031 

-7 

.9273 

.8U57 

1.U637 

5.7615 

.038 

-8 

1.0520 

1.1072 

2.U988 

12.2010 

.033 

-9 

1.0323 

1.2839 

3.8978 

23.652U 

.016 

-10 

.9751 

.8368 

1.3691 

5.3359 

.022         2270  microsecs 

Theor. 

1.0 

1.0 

2.0 

9.0 

.039  (1<#) 

-11 

.985U 

1.0562 

2.U39U 

11.612U 

.018 

-12 

.9998 

1.0U95 

2.2U30 

10.3795 

.OlU 

-13 

.98U1 

1.0UU7 

2.2683 

10.1315 

.037 

-111 

1.020U 

1,0517 

1.9077 

7.1153 

.021 

-15 

.9513 

.7981 

1.186U 

U.1650 

.022 

-16 

.9592 

.9722 

1.7937 

6.9373 

.oia 

-17 

1.0U03 

I.IO69 

2.2U5U 

9.935U 

.017 

-18 

1.0329 

1.0635 

1.9U97 

7.6527 

.018 

-19 

.9711 

.901U 

1.6697 

6.832U 

.022 

-20 

.9835 

.9U51 

1.7Uli5 

7.1993 

.020 

-21 

1.0051 

1.0281 

2.15U7 

9.U030 

.015 

-22 

.9631 

.8759 

1.3956 

U.915U 

.020 

-23 

.9951 

.9813 

1.778U 

7.1579 

.015 

-2U 

1.0176 

1.1237 

2.5168 

12 .1998 

.020 

-25 

1.0278 

1.066U 

2.2520 

10.2068 

.025 

-26 

1.0317 

1.0396 

2.1303 

9.3136 

.027 

-27 

1.0188 

.9595 

1.6659 

6.7573 

.02U 

-28 

.933U 

.8338 

1.3261 

U.7659 

.038 

-29 

.9U68 

.9008 

1.8U11 

8.U889 

.03U 

-30 

1.0172 

.9759 

I.80I4I 

7.57U6 

.029 
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Addendum  1„ 

Results  of  the  investigation  of  the  tails  of  the  Marsaglia  normal 

random  number  generator. 

The  test  results  belie  the  idea  resulting  from  investigations 
made  in  Section  U.3.  A  typical  series  of  results  of  the  tests  are 
shown  below.  The  first  line  gives  the  theoretical  cumulative  sample 
result  based  on  a  sample  size  of  100,  The  following  data  shows 
the  experimental  results.  The  first  interval  is  from  minus  infinity 
to  -2,56,  following  intervals  are  0,16  in  width.  Thus  only  the 
negative  half  of  the  distribution  is  shown. 
.52  .82  1.25  1.88  2.7U  3.92  5.U8  7.U9  10.03  13. lU  16.85  21.19  26.11 
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